^ Received: may 24, 1999, Accepted: August 19, 1999 

HYPER VERSION 



Implications of a new solar system population 
of neutralinos on indirect detection rates 



Lars Bergstr6m a , Thibault Damour 6 , Joakim Edsj6 a , 
Lawrence M. Krauss c and Piero Ullio a 

a Dept. of Physics, Stockholm University, Box 6730, SE-1 13 85 Stockholm, Sweden 
\ b Institut des Hautes Etudes Scientifiques, 

00 



EE 

35 route de Chartres 91440 Bures sur Yvette, France 

Departments of Physics and Astronomy, Case Western Reserve University 
10900 Euclid Ave, Cleveland OH 44106-7079 

in . 
o. 

Abstract: Recently, a new Solar System population of weakly interacting massive 
. particle (WIMP) dark matter has been proposed to exist. We investigate the impli- 

cations of this population on indirect signals in neutrino telescopes (due to WIMP 
Oh. annihilations in the Earth) for the case when the WIMP is the lightest neutralino 

of the MSSM, the minimal supersymmetric extension of the standard model. The 
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velocity distribution and capture rate of this new population is evaluated and the 
flux of neutrino-induced muons from the center of the Earth in neutrino telescopes 
is calculated. The strength of the signal is very sensitive to the velocity distribution 
of the new population. We analytically estimate this distribution using the approx- 
imate conservation of the component of the WIMP angular momentum orthogonal 
to the ecliptic plane. The non-linear problem of combining a fixed capture rate 
from the standard galactic WIMP population with one rising linearly with time from 
the new population to obtain the present-day annihilation rate in the Earth is also 
solved analytically. We show that the effects of the new population can be crucial for 
masses below around 150 GeV, where enhancements of the predicted muon flux from 
the center of the Earth by up to a factor of 100 compared to previously published 
estimates occur. As a result of the new WIMP population, the next generation of 
neutrino telescopes should be able to probe a much larger region of parameter space 
in the mass range 60-130 GeV. 
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1. Introduction 

It has been known for almost 15 years that Weakly Interacting Massive Particles 
(WIMPs) can elastically scatter inside the Sun and Earth, leading to their subsequent 
capture and annihilation in the cores of these objects, producing an indirect neutrino 
signature that might be accessible to neutrino telescopes |]. Recently, however, it 
has been demonstrated that the scattering process in the Sun can populate orbits 
which subsequently result in a bound Solar System population of WIMPs p], || and 
which can be comparable in spectral density, in the region of the Earth, to the 
Galactic halo WIMP population. This new population consists of WIMPs that have 
scattered in the outer layers of the Sun and due to perturbations by the other planets 
(mainly Jupiter) evolve into bound orbits which do not cross the Sun but do cross 
the Earth's orbit. This population of WIMPs should have a completely different 
velocity distribution than halo WIMPs and will thus have quite different capture 
probabilities in the Earth. The predicted WIMP abundance, and spectrum, relevant 
for direct detection have been calculated [0, [3j. Here we focus on capture in the 
Earth, and the predicted indirect neutrino signature. Other studies of solar system 
populations of WIMPs can be found in |], |], 

Following [0, |[, we consider here a special WIMP candidate, the neutralino, 
which arises naturally in supersymmetric extensions of the standard model. We 
evaluate the capture of neutralinos and the resulting neutrino-induced muon flux 
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within the Minimal Supersymmetric extension of the Standard Model (MSSM) (for 
a review of neutralino dark matter, see ref. We find that this flux is sensitively 
dependent upon the Solar System WIMP velocity distribution, and can exceed by 
two orders of magnitude that predicted for Galactic halo WIMPs. Although our 
numerical results apply only to supersymmetric WIMPs, we expect our qualitative 
conclusions to be generally valid for any WIMP model of dark matter. 



2. Definition of the supersymmetric model 

We work in the Minimal Supersymmetric Standard Model (MSSM). In general, the 
MSSM has many free parameters, but with customary assumptions || (which do 
not significantly affect the general results of this paper) we can reduce the number of 
parameters to seven: the Higgsino mass parameter /x, the gaugino mass parameter 
M 2 , the ratio of the Higgs vacuum expectation values tan (3, the mass of the CP-odd 
Higgs boson m^, the scalar mass parameter mo and the trilinear soft SUSY-breaking 
parameters A& and A t for the third generation. In particular, we do not impose any 
restrictions from supergravity other than gaugino mass unification. We have made 
some scans in parameter space without the GUT relation 3Mi = 5M 2 tan 2 9w for 
the gaugino mass parameters Mi and M 2 . This mainly has the effect of allowing 
lower neutralino masses to escape the LEP bounds and is not very important for 
this study. Hence, the GUT relation for M 1 and M 2 is kept throughout this paper. 
For a more detailed definition of the parameters and a full set of Feynman rules, 
see refs. [7L EJ. Of course we assume i?-parity conservation, which ensures that the 
lightest supersymmetric particle is stable and therefore an excellent dark matter 
candidate. Its relic abundance today is fixed by its freeze-out abundance in the early 
Universe, which in turn is given by the interplay between the interaction rate and 
the expansion rate at a temperature around T ~ m x /20, where m x is the neutralino 
mass (i.e. it is non-relativistic at freeze-out and would act as cold dark matter in 
structure formation). 

The lightest (and, with i?-parity conservation, stable) supersymmetric particle is 
in most models the lightest neutralino, which is a superposition of the superpartners 
of the gauge and Higgs fields, 

X? = N n B + N 12 W 3 + N 13 H° + N U H° 2 . (2.1) 

It is convenient to define the gaugino fraction of the lightest neutralino, 

Z g =\N u \ 2 + \N 12 \ 2 . (2.2) 

For the masses of the neutralinos and charginos we use the one-loop corrections 
as given in ref. and for the Higgs boson masses we use the leading log two- 
loop radiative corrections, calculated within the effective potential approach given in 
ref. m. 



Parameter 




M 2 


tan/? 


m A 


m 


A b /m 


A t /m 


Unit 


GeV 


GeV 


1 


GeV 


GeV 


1 


1 


Min 


-50000 


-50000 


1.0 





100 


-3 


-3 


Max 


50000 


50000 


60.0 


10000 


30000 


3 


3 



Table 1: The ranges of parameter values used in our scans of the MSSM parameter space. 
Note that several special scans aimed at interesting regions of the parameter space have 
been performed. In total we have generated approximately 1.4 x 10 5 models which are not 
excluded by accelerator searches. 



We have made extensive scans of the model parameter space, some general and 
some specialized to interesting regions, where high rates are possible. In total we 
have made 33 different scans of the parameter space. The scans were done randomly 
and were mostly distributed logarithmically in the mass parameters and in tan/3. 
For some scans the logarithmic scan in \x and M 2 was replaced by a logarithmic scan 
in the more physical parameters m x and Z g /{1 — Z g ), where m x is the neutralino 
mass. Combining all the scans, the overall limits of the seven MSSM parameters we 
use are given in table [[]. 

We check each model to see if it is excluded by the most recent accelerator 
constraints, of which the most important ones are the LEP bounds [IT] on the lightest 
chargino mass, 

[91 GeV, |m + — m Y o \ > 4 GeV , , 

m + > < Xl 1 (2 3) 

Xl [ 85 GeV , otherwise v ' ; 

and on the lightest Higgs boson mass m H o (which range from 72.2-88.0 GeV de- 
pending on sin(/3 — a) with a being the Higgs mixing angle) and the constrains from 



b — > S7 [JT2 ] . For each allowed model we calculate the relic density of neutralinos 



Q x h 2 , where Q x is the density in units of the critical density and h is the present 
Hubble constant in units of 100 km s _1 Mpc -1 . We use the formalism of ref. |L3 



for resonant annihilations, threshold effects, and finite widths of unstable particles 
and we include all two-body tree-level annihilation channels of neutralinos. We also 
include so-called coannihilation processes in the relic density calculation according 
to the analysis of Edsjo and Gondolo M. 

Present observations favor h = 0.65 ±0.1, and a total matter density Qm = 
0.3 ±0.1, of which baryons may contribute 0.02 to 0.08 Not to be overly 

restrictive, we accept Q x h 2 in the range from 0.025 to 1 as cosmologically interesting. 
The lower bound is somewhat arbitrary as there may be several different components 
of non-baryonic dark matter, but we demand that neutralinos are at least as abundant 
as required to make up the dark halos of galaxies. In principle, neutralinos with 
Q x h 2 < 0.025 would still be relic particles, but only making up a small fraction of 
the dark matter of the Universe, so we do not consider this alternative here. 



3. Speed distribution of the new neutralino population 

The distribution of Solar System WIMPs relevant to terrestrial capture arises from 
scattering of Galactic halo WIMPs into highly elliptical Solar System orbits, with 
semi-major axes in the range a ea rth < a < 2.6a eart h- These are then perturbed by 
Jupiter and the other planets into orbits which do not intersect the Sun. Analytical 
methods were used in refs. |2|, [| to estimate the basic features of this distribution. 
Here we shall relax an approximation (exactly radial orbits) made in refs. |2], |3| and 
derive an improved analytical estimate of this distribution. We nevertheless expect, 
because of their origin in Solar scattering, that most of these orbits will be close to 
radial in the vicinity of the Earth. 

In ref . || , the phase space distribution of WIMPs at the Earth was investigated 
including the effects of the Earth being deep in the Sun's potential well. There it was 
found that even though the velocity distribution of WIMPs at the Earth is different 
than it would be in free space, both Jupiter, the Earth and Venus will disturb the 
orbits of these unbound WIMPs into bound orbits that have the same phase space 
distribution as would be the case in free space. This relies on the fact that the time- 
scales for transferring WIMPs between these unbound and bound orbits are shorter 
than the age of the Earth so that an equilibrium occurs. This is valid for low velocities 
relative to the Earth only (which are the most important ones for capture) and it was 
found that for some higher-velocity orbits (velocities > 27 km/s), the time-scale for 
filling or depleting these bound orbits is too long. It turns out that the (nearly) radial 
orbits of interest for this work happen to be in precisely this regime. More precisely, 
our new population is distributed along a thin vertical "needle" , located within the 
"unfilled bound orbits" region of fig. 3 of ref. |J, starting at u = 30 km/s on the 
left part of the horizontal axis. This is of course an advantage for detection since it 
means that these WIMPs will stay in these orbits without being much affected by 
diffusion within the Solar System. 

The detection rates we discuss below are very sensitive to the actual form of 
the velocity distribution, especially at low velocities (with respect to the Earth). 
The reason for this is easy to understand from simple kinematics. A very heavy 
WIMP (heavier than iron) which encounters the Earth cannot be stopped by a single 
encounter if its velocity is large. (And the optical depth for repeated scatterings is 
very small.) Therefore, it is the strength of the distribution at low energies which 
determines the capture rate and therefore the muon neutrino rate from the center 
of the Earth. An exception to the rule of only low velocities being important occurs 
if the WIMP has a mass which closely matches one of the most abundant elements 
of the Earth. Then kinematics allows all the kinetic energy of the WIMP to be 
transferred to a nucleus in a single collision. As we shall discuss in detail below, 
conservation of the z component of the angular momentum constrains very much the 
component of the velocity of the WIMPs along the motion of the Earth, and thereby 



the lowest attainable speed with respect to the Earth. As a consequence, we shall 
see that only neutralinos of the new population lighter than around 150 GeV will be 
captured by the Earth. 

Let us now turn to the best estimate we can presently make of the velocity 
distribution, as seen in an Earth frame, of the new neutralino population. In refs. 0, 
|, as a first approximation, a model was used where the WIMP orbits were nearly 
radial (i.e. with a fixed eccentricity, very near e = 1). Here, we shall go beyond this 
approximation by taking into account the distribution in the z component of the 
angular momentum (which will turn out to be the crucial one for our purpose). Let 
us first define our notation. 

In an heliocentric frame, the WIMP phase-space distribution function would 
read dN = /(x, v) <i 3 x A <i 3 v, where x denotes the WIMP heliocentric position and 
v its heliocentric velocity. We could start from /(x, v), written as a function of the 
integrals of motion of the WIMP (so that it satisfies Liouville's equation), to derive 
the capture rate in the Earth. It is, however, more convenient to take advantage 
from the start of the simplification brought by the spherical symmetry of the Earth. 
Let u = v — V£ denote the incoming velocity (far from the Earth) of the WIMP in 
an Earth reference frame (here denotes the heliocentric velocity of the Earth). 
The WIMP capture rate by the spherically symmetric Earth depends only on the 
angular average of the velocity distribution, and therefore is a function only of the 
distribution of "speeds" u = |u|. We can write the speed distribution of the local 
WIMP number density as 

dn(u) = n(ai) dn(u) , (3.1) 

where n(ai), with a\ = 1AU denoting the radius of the Earth's orbit, is the total 
space density around the Earth (WIMPs per cm 3 ), and where dn(u) is the normalized 
speed distribution: J dn(u) = 1. The integrated WIMP space density ra(ai) has been 
estimated in refs. [0, [|. We shall quote the result below when we need it. We focus 
first on the normalized speed distribution dn{u). 

At any moment we can introduce an heliocentric spatial reference frame (x, y, z) 
such that the x-axis is in the Sun-Earth radial direction, and the |/-axis is directed 
along the orbital motion of the Earth. In this frame 

u 2 = {y-w E f = w 2 + v 2 E -2v E v y . (3.2) 

As we are considering a fixed value r = |x| = a\ of the WIMP radial position, the 
squared heliocentric velocity v 2 can be expressed in terms of the semi-major axis of 
the WIMP orbit: 

v 2 (r = ai ) = v\ (2 - ^) . (3.3) 

On the other hand, the longitudinal ( "along the track" ) velocity v y entering eq. ( |3.2j ) 
can be expressed in terms of the ^-component J z of the angular momentum of the 
WIMP, v y = J z /a\ (note that the z axis is orthogonal to the ecliptic plane). Finally, 



the (square of the) local speed of the WIMP can be entirely expressed in terms 
of the two adiabatic invariants a and J z (i.e. the two Delaunay action variables 
L = \/Gn M Sun a and H = yjG^ M Sun a(l — e 2 ) cosz in the notation of ref. 0): 



2 2/o "1 

u — v E I 3 



2J Z 



(3.4) 



ai v E 



We recall that ve = (Gn ^sun/^i) 1 ^ 2 = 27.98 km s _1 denotes the Earth orbital 
speed, and ai = r E = 1 AU. 

Let us recall from ref. || that, to a good approximation, the two WIMP action 
variables L and H are secularly conserved (modulo some small random- walk-type 
diffusion caused by near encounters with the inner planets), while the third action 
variable G = \J Gn Ms un a(l — e 2 ) will secularly evolve under the influence of plan- 
etary perturbations. [The averaged planetary Hamiltonian 7i p (Kozai Hamiltonian) 
was discussed in section 3 of 0.] The evolution of G means that the eccentricity of 
the WIMP orbit oscillates between values very near one and smaller values. We see 
from eq. Q3.4p that these eccentricity oscillations have no effect on the local speed 
distribution. The latter depends only on the well conserved quantities L and H, or 
a and J z . We can then estimate the present speed distribution by means of the ini- 
tial distribution of a and J z . ["Initial" means here at the time when the considered 
WIMP of the new population exited from the Sun after having been captured by 
scattering on a nucleus in the outskirts of the Sun.] The initial ^-component of the 
WIMP angular momentum was 



Here, ao(l — eo) is the initial perihelion distance which is approximately equal to 
the (effective) Sun radius Rs = 0.907 -Rsun- The remaining factor 1 + eo can be 
approximated by 2, so that 



where the cosine of the initial inclination can be (approximately) considered as being 
a random variable uniformly distributed between —1 and +1. 

Because of the overall (approximate) axial symmetry of the Solar System, we 
expect the conservation law of J z to be rather accurate: J^ ow ~ J°. However, as 
the nonconservation of J z would have a crucial effect on the capture rate of WIMPs 
by the Earth (by allowing for lower speeds in an Earth frame, through the last term 
in eq. ( J3.4Q ) we shall introduce a phenomeno logical parameter A to parametrize a 
possible nonconservation of J z . More precisely, we shall assume that 



J z = ^/G N M Sun a (l - e )(l + e ) cosi • 



(3.5) 




(3.6) 



now 



A J? 



(3.7) 
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and study not only the "standard" case where A = 1, but also the case where, say, 
A = 2, which corresponds to a 100% violation of the conservation of J z . Then, 
eq. ( |3.4| ) reads 

u 2 = v% (3 — e cosio) , (3.8) 

where 



2 R 

2X\I -~ 0.18377 A. (3.9) 

Eq. (|3.8|) tells us that the distribution function of the squared speed u 2 is obtained 
by taking the convolution between the distribution of the variable b = a\ja (which 
is essentially the WIMP heliocentric energy) and the (uniform) distribution of the 
variable c = e cosi - If we denote b' = b + e cosz , the distribution ip'(b') db' of b' is 
given in terms of the distribution ip(b) db of b = a%/a by 

The normalized distribution (p(b) db (with J tp(b) db = 1) has been derived in 
ref. 0. It is written in eq. (6.7) there, in terms of the related variable x = 2a/a\ = 
2/6. (Ref. U made the approximation of nearly radial orbits in deriving eq. (6.7), but 
this approximation is not crucial for our present purpose. The crucial new feature 
for the present work is the allowance for the along-the-track term proportional to J z 



in eq. (|3.4 ) which is the only one which can significantly affect neutralino capture by 



the Earth.) Transcribing eq. (6.7) of in terms of the variable b = 2/x yields 

V( b ) db = N » fro.a { 2 b _ 6) i/2 9 ( b - 2^) *(2 " 6) . (3-11) 
with the normalization constant 



1O.I 



N b = — = 0.456579 . (3.12) 

4(5.2) 1 ; 

Here, as above, 6{x) denotes the step function (9(x) = 1 if x > 0, 9(x) = if x < 0). 

Eq. ( |3.10| ) then gives the distribution in b', or equivalently in the scaled squared 
speed 

^2 ( u 



U = u 2 =y—j =3-b' = 3-b-e cosi . (3.13) 

Before the smearing of eq. (|3.10|) (due to the convolution with the J z -distribution) 
the original distribution (p(b) db of eq. ( p.ll ) featured both a cut-off at b\ — a\/a — 2 



and a cut-off at 62 = 1/2.6. The former cut-off corresponds to a lower- velocity cut-off 
at u — Ve(3 — bi) 1 / 2 = Ve = 29.78 km s _1 . (It corresponds to radial orbits which 
barely reach up to the Earth orbit.) In our present refined treatment, where we take 



into account the along-the-track motion linked to the distribution in J z , this lower- 
velocity cut-off will be shifted to a lower value (coming from a positive J z in eq. ( p.4|) ). 
For capture by the Earth (which has a small escape velocity, ranging between about 
11 and 15 km s -1 ) the low-velocity part of the speed distribution is the only one 
to play a crucial role. We can take advantage of this fact to derive a simplified 
analytical expression of the speed distribution dn(U) = <p'(3 — U) dU = <p'(b')db'. 
Indeed, we can approximate the smearing, eq. Q3.1UI ), as acting only on the crucial 
low-velocity cut-off factor 9(2 — b)/(2 — b) 1 ' 2 in eq. fl3.11|) . The smearing of the 



smoothly varying factor l/b os would only introduce a small fractional correction of 
order e 2 ~ 3%. With this approximation the smearing integral of eq. ( p.lOj ) can be 
explicitly computed. The final result for the speed distribution reads 

dn = N'J £ (U)dU = N' £ f £ (u 2 )2udu, (3.14) 

where we recall that U = u 2 = (u/ve) 2 , where N' e is a normalization constant (equal 
to Nf,, eq. ([3.12 ), modulo 0(e 2 ) fractional corrections 1 ) and where f E {U) is explicitly 



given by 

f^ U ) = H ^ X u)o U J - £ [VU-l + ee(U-l+e)-VU-l-e6(U-l-s 

(3.15) 

Here f/ max = 3 — 1/2.6 ~ 2.6154 is the upper cut-off, and e is given by eq. ( |3.9D (with 
A = 1 being our standard estimate, and A = 2 being an extreme phenomenological 
possibility). Note that the speed distribution extends only between a minimum cut- 
off at 

Mmin = v E Vl - e (3.16) 

and a maximum at u max = ve V^max — 48.16 km s _1 . (The exact maximum cut-off 
given by eq. ( |3.10| ) is modified by 0(e) terms, but this modification is unimportant 
for our present purpose.) The numerical value of the lower speed cut-off, eq. (|3.16| ), 
is 26.905 km s -1 for the standard case A = 1, and 23.683 km s _1 for the extreme 
case A = 2. The precise value of this lower speed cut-off is important because it 
translates into an upper cut-off in the masses of the WIMPs that can be captured by 
the Earth. Indeed, as will be recalled below, the kinematics of the maximum energy 
loss in the collision of a WIMP of mass mx with a nucleus of mass rriA located at 
radius r in the Earth (with local escape velocity f eS c(^)) is such that capture is only 
possible if 

(m x -m A y \v esc (r)J 



1 In our calculations we used N' e = N b = 0.456579. We verified numerically that this induces an 
error of less than 2% for A = 1 . 
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This implies the following upper limit on the mass of a WIMP capturable by the 
Earth 

m x < (l + 2w + 2 Vw 2 +wjm A , (3.18) 

where we denoted w = (v esc (r) / u m i n ) 2 . The highest capturable mass is reached when 
the WIMP scatters on an iron nucleus (m^ = mp e — 52 GeV) located at the center 
of the Earth (where t> eS c(0) — 15 km s _1 ). In our standard case (A = 1) where 
u min ~ 26.905 km s -1 this yields the upper limit 

m x < 2.90 m Fc ~ 150 GeV, (3.19) 

on the other hand in the case A = 2 the limit is shifted to about 170 GeV. 



4. Overdensity of the new population 

We follow the notation of 0, |3| and write the contribution from this new population 
of WIMPs to the usual halo WIMP density as 



n(a>i) _ (secondary) WIMP density at the Earth 
rix halo WIMP density at infinity 



(4.1) 



where 



5.44 x 10 36 _ „ _ 2 0.212 ( _ 10) (A . 

5 E = r- x 9tot GeV cm 2 = ^ gL t ' . (4.2) 

K/220kms- 1 ) y (V220kms J 

Here, g^ t 10) = 10 10 5f to t(GeV) 3 , and g tot = J2a (/a/ m i) °a K s a , where f A is the mass 
fraction of element A in the Sun, and K S A is the surface value of the capture function 
on the element of mass number A in the Sun according to eq. (2.20) in ref. 0. 

The scattering rate of WIMPs in the outer layers of the Sun (which causes the 
fast halo WIMPs to lose enough energy to enter bound orbits close to the Earth's 
orbit) is proportional to <JaK a , which can be calculated once the parameters of the 
SUSY WIMP in question are fixed. (For the elemental abundances in the Sun, we 
use the compilation of Bahcall and Pinsonneault ||15|| .) 

In fig. m we display the values of g^ 10 ' 1 versus neutralino mass for our set of 
supersymmetric models. As can be seen, in some cases values approaching or even 
exceeding unity can be obtained (confirming the results of 0, ||). The spread is 
very large, however, and some models give orders of magnitude smaller values. As 
would be expected, the models with the highest values of g[o^ are the same models 
which give high scattering rates in direct detection experiments (the value of the 
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Figure 1 : The value of the parameter g[ ot 10 ^ , which is related to the scattering probability 
near the Solar surface, as a function of neutralino mass. The different symbols represent 
different ranges for the spin-independent cross section, where o^™ represents the most 
stringent claimed upper limit on this cross section for an assumed halo density of 0.3 
GeV/cm 3 , given in [16] (see text for a discussion). 



local halo density is assumed to be 0.3 GeV cm -3 ). This can be seen by the coding 
of the displayed points in terms of one claimed limit ||16| on the spin-independent 
scattering cross section 2 . 

Note that the highest values of g^ot 10 "* allowed by this claimed limit on the scat- 
tering rates is g^/max — 0.4. Such values imply only a rather modest (but not 
negligible, and important as independent check) effect of the new population in di- 
rect experiments (see refs. @, U). On the other hand, we shall see below that they 
can imply a large effect on indirect detection rates. 



2 Note that in ref. ||, || larger values of <?tot were obtained because the possibility of smaller halo 
WIMP densities was allowed, and also because more conservative claimed limits were used. Here 
instead, in order to be conservative in our estimate of predicted muon fluxes in detectors, we use 
as our fiducial value of the WIMP scattering cross section the most stringent claimed cross section 



upper limit 1 16 



10 



5. Capture rate of WIMPs by the Earth 



The total space density of WIMPs, per speed element, in the vicinity of the Earth 
reads 



dn°x + dn 



new 
X 



nx dn old (u) + ra new (ai) dn new (u) 
n x [dn old (u) + 5 E dn new (u)] . 



(5.1) 



The "old" distribution dn old (u) can, to a good approximation ||, be considered to 
be simply the incoming Galactic one, shifted to the frame of the moving Earth. 
The "new" (secondary) distribution dn new (u) is that defined by eq. (|3.14|) above. 
The present overdensity ratio 5e is that discussed in section f| above. Let us now 
write down the capture rate of WIMPs by the Earth |17| , |18| . We shall utilize the 
result of Gould 



18 1, using the notation of ref. 



As said above, we work here 



with the angle-averaged speed distribution, which reads, in the notation of ref. 
dnx{u) = in u 2 du f ^(w) . Using eqs. (2.20)-(2.22) of ref . 0, considered for J min - 
we get for the capture rate on element A: 



0. 



1 / \ dnxiu) 

dN A = -j d 3 xn A (r)a A 



u 



exp 



(u + a) 



2Q; 



9 a da , 



with 



Here 



9 n = 6 



(r)- — 



4mx rriA 



u 



^ (m x ± m A ) 2 ' 



(5.2) 



(5.3) 



(5.4) 



a A would be the total scattering cross section on nucleus m A in absence of the form 
factor F A (Q) = exp (—Q/Q A ), with Q = E heiore — E after denoting the energy loss and 
Q A the standard value of the coherence energy 0, recalled in eq. (2.28) of 0, and 
a is linked to the semi-major axis of the WIMP, after capture by the Earth, by 



a = 



Gn ^Earth 



(5.5) 



Note that we work here in a geocentric frame (with u = v^? — v^). The differential 
capture rate, eq. (|5.2|) , must be integrated over a and the volume of the Earth, and 
must be summed over the label A denoting the various elements in the Earth. The 
integral over a is limited by the theta function 8 a , eq. (|5.3| ), and by the fact that 
a must be positive (which expresses the fact that the WIMP lost enough energy by 
scattering on element A to end up being bound to the Earth). These two constraints 
restrict the range of variation of a to < a < a m (r,u), with 



a. 



,(r,u) =(3*v 2 esc (r) 



0t 



u 



(5.6) 
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Region 


fo f~Mg fsi /Fc 


r < 3488 km 

3488 km< r < 6378 km 


0.05 0.90 
0.44 0.23 0.21 0.059 



Table 2: The mass fractions, f+ of various elements in the Earth. Data from ref. [19]. 



and imply the presence of a remaining theta function restricting the allowed values 
of u: 

O(^vUr)-u'). (5.7) 



On, 



Finally, integration over a (followed by integration over <i 3 x and drix{u), and sum- 
mation over A) leads to a total capture rate 



E 2Qa 



,4 
X 



(3+mx 
[exp (- 



d 3 xn A (r) a A 



dnx{u) 



x 



u 



(5. 



This capture rate is a linear func- 
tional of the WIMP spectrum dnx{u). 
Therefore the decomposition, eq. (|5.1|) , 
implies a corresponding linear decom- 
position of the capture rate 



C x 



C o\d + (J* 




(5.9) 

To evaluate C new we have to perform 
the integral in eq. Q5.8p explicitly. Avail- 
able analytical results in the literatu- 
re do not apply since they assume a 
Maxwell-Boltzmann velocity distribu- 
tion. Since we have spherical symme- 
try, we are left with a double integral, 
one over the distance from the Earth's 
center, r, and one over the velocity u. 
For the composition and density of the 

Earth as a function of r we use the distributions given in ref. [|BJ . From the density 
distribution, we can calculate the gravitational potential as a function of r and hence 
the escape velocity v csc {r). The density and escape velocity as a function of r are 
given in fig. |2|. In table ^] we show the mass fraction of the (for our purpose) most 
important elements for different regions in the Earth. 



4000 5000 6000 
Radius (km) 



Figure 2: Density (a), and escape velocity 
(b), as functions of the distance from the cen- 
ter of the Earth. 
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6. Neutralino annihilation rate in the Earth 



Let N(t) be the total number of neutralinos trapped, at time t, in the core of the 
Earth. The annihilation rate of neutralino pairs can be written as 

T a {t)= l -C a N\t). (6.1) 

The evolution of N(t) is the result of the competition between capture and 
annihilation: 

^- = C tot (t)-C a N 2 (6.2) 

(the factor 2 difference in the annihilation terms in eqs. (|6.1|) and fl6.2Q comes from the 
fact that neutralinos annihilate in pairs). The constant C a entering equations (|6.1 ) 
and ( |6.2| ) is linked to the annihilation cross- section a a , and to some effective volumes 
Vj, j = 1, 2, taking into account the quasi-thermal distribution of neutralinos in the 
Earth core: 

V 2 

C a =(a a v)y^, (6.3) 

/ ■ \ -3/2 

2 ' 3 x 1025 [w^vj cm3 - (6 ' 4) 

Usually, the capture rate C in eq. ( |6.2| ) is time- independent, being given by the "old", 
Galactic WIMP population: C usual = C old = const. In that case, eq. O can be 
simply solved by separating the two variables N and t, i.e. 

' = }* = } oh=c3H- (6 ' 5) 

This yields for the standard, "old" annihilation rate 

rf = \c a N 2 = ^ old tanh 2 ( 7 1) , (6.6) 

with 

7 ^(C old C a ) 1 /2. (6 . 7 ) 

In the case considered here, where capture comes both from the direct Galactic 
population and from the secondary population discussed in refs. @, |J, the total 
capture rate increases linearly with time, because of the linear-in-time build up of 
the new population: 

C tot (t) =C old + Ct, (6.8) 

where 

/'"mew 

C=-—. (6.9) 

Here, C new is the present value of the capture rate of the new population (that 
discussed in the previous section which assumed an overdensity Se built up over the 
full time ts) and ts = 4.5 Gyr is the build up time, i.e. the lifetime of the Solar 
System. 



13 



Note that the problem of determining the total, present annihilation rate, in 
presence of the new population, is a non-linear one. The answer cannot be simply 
split as the effect of the old population, plus the effect of the new one. The differen- 
tial equation to be solved, eq. ( |6.2|) with eq. (|6~8|), is a Riccati equation. In general, 
such an equation cannot be solved analytically. However, in the present case (cor- 
responding to the restricted type of equations originally discussed by Riccati), the 
problem can be solved in terms of Airy functions. To do that let us first scale the 
variables to simplify the Riccati equation. Let us introduce the reduced variables x 
and y by posing 

C tot (t) = C oXd + Ct = ax, N = 0y, (6.10) 

where 

a = C 2 ' 3 C7 1 / 3 , = C 1 / 3 C- 2 ' 3 (6.11) 

(the notation a here should not be confused with our unrelated previous use of this 
letter). In terms of x and y the evolution equation ( |6.2|) reads simply 

^- + y 2 = x. (6.12) 

ax 

If we now introduce a new variable u(x) by posing y(x) = (du/dx)/u (this use 
of the letter u should not be confused with the notation above for the speed with 
respect to the Earth), eq. (|6.12|) becomes the Airy equation: 

^=xu. (6.13) 

The general solution of ( |6.13|) reads u(x) = aAi(a;) + bBi(x), where Ai (x) is the 
usual Airy function (decreasing for x —>■ +00) and Bi (x) the complementary solution 
(increasing as x —>■ +00). The boundary condition of our problem is that y vanishes 
when x takes it initial (positive) value 

/-<old 

Xi = (6.14) 

a 

(corresponding to t — 0). We are then interested in the final value yj = y(xf) of y 
corresponding to the final (present) value of x: 

C tot (ts) C old + Ct s 
xt = = . 6.15 

a a 

This final value is given by 

Ai'^Bi'^-Bi'Qr^Ai'^) 
Vf Ai(x 7 ) Bi'fa) - Bi(x f ) Ai'(xi) ' l ' ' 

In terms of yj the total, present annihilation rate reads 

K ot (t s ) = lay 2 f . (6.17) 
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Figure 3: The enhancement factor, eq. ( 6.1S| ), as a function of the neutralino mass. The 
points are symbol coded according to whether they exceed the claimed direct detection 
bound a si [0] from the Dama experiment (squares), give direct detection rates within a 
factor of 10 of those limits (pluses) or give rates more than 10 times smaller than these 
limits (crosses). 



It is useful to consider also the amplification factor A brought by the new pop- 
ulation, i.e. the ratio 

r tot (*<?) ay 2 f 

A = 1 + S = ° \ A S) = — -, (6.18 

T° ld C° ld tanh 2 ( 7 t s ) 

where we have also denned the enhancement factor £ = A — 1. Note that since 
the muon flux is directly proportional to the annihilation rate, $ M = k T a , with k 
depending on the neutralino mass and the 'hardness' of the annihilation spectrum, 
the enhancement factor reads 

-ptot _ -pold (jjold+ncw _ ^old 

£ = ^Y^= ~ $old ~ • (6-19) 

^ a fj. 

In fig. [| the enhancement factor for our set of MSSM models is shown against 
the neutralino mass. As anticipated, the amplification mechanism is only operative 
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for masses smaller than around 150 GeV. Also, for small masses (less than around 
60 GeV) the enhancement is relatively small. This is due to a combination of fac- 
tors, one being that an important contributor to the scattering in the Sun which 
gives rise to the radial, "planetary" orbits is iron (in fact, it is important also for 
the capture in the Earth), and for WIMP masses smaller than iron the momen- 
tum transfer is less effective. In the WIMP mass region between 60 and 130 GeV, 
on the other hand, enhancements up to a factor of 100 are possible. This can po- 
tentially be very important for the discovery potential for dark matter of neutrino 
telescopes. 

The points in the figure are symbol-coded according to what direct detection 
rates they correspond to. In particular, the squares indicate neutralino models which 



exceed the claimed Dama experiment upper limit |16| , assuming a halo WIMP den- 
sity of 0.3 GeV/ cm 3 . In computing the limits from Dama we have included the effects 
of the new population, but because of the low velocities their effect on present-day 
detectors with their relatively high recoil energy thresholds is quite small, and in fact 
hardly visible in our figures. 

7. Muon fluxes 

We calculate the resulting neutrino-induced muon fluxes essentially as described 



in 



20]. The decay and/or hadronization of the annihilation products as well as the 
neutrino interactions are simulated with Pythia 6.115 [^TJ. All two-body annihi- 
lation final states (at tree level) are included as well as the 1-loop induced 2 gluon 
and Z gluon. We use a value of 0.3 GeV/cm 3 for the local neutralino mass density 
(leaving a study of the change in our analysis caused by variations in this and other 
astrophysical parameters to future work). 

We compare the predicted rates with the present experimental upper bound, 



where we have taken the Baksan results |E§ as a representative (the MACRO p3 



limits are very similar). The muon threshold is set to 1 GeV. In fig. |] we show our 
predicted total muon rates from the direction of the center of the Earth, compared 
to the rates from the "old" population alone. (Note, that due to the non-linearity of 
the rate equations, one cannot just sum the contributions of the "new" flux and the 
"old" flux separately.) As can be seen, enhancements between one and two orders 
of magnitude are frequent (even though the enhancement for the models with the 
highest predicted absolute flux is generally somewhat smaller). 

In fig. § we show the predicted absolute fluxes against neutralino mass for our 
sample of models, without and with the new population. The Baksan limit is shown 
as the nearly horizontal line. As can be seen, the effect on this scatter plot of the 
new population is striking. Many models around 60-130 GeV have been boosted 
up to higher fluxes. Several are now above the Baksan limit and others are within 
detectability in the near future. Around 80 GeV, the increase in flux for models 
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Figure 4: The muon flux including both the old and the new population versus the flux 
obtained if only the old contribution is included. The points are symbol coded according to 
whether they exceed the claimed direct detection bound <jsi [0] from the Dama experiment 
(squares), give direct detection rates within a factor of 10 of those limits (pluses) or give 
rates more than 10 times smaller than these limits (crosses). 



with high fluxes can be as high as a factor of 10 (for low flux models, the in- 
crease can be even higher as seen already in fig. One may note that some of the 
models previously thought to be allowed by Baksan but above the claimed Dama 
bound, should really be considered as ruled out by Baksan with even greater confi- 
dence. 

We have investigated the composition of the neutralinos that cause the largest 
enhancement of the muon flux and that at the same time have high absolute fluxes. 
We find that they are generally of mixed type (neither very pure Higgsino nor gaug- 
ino), and that the biggest enhancements occur for neutralinos which simultaneously 
have large spin-dependent and spin-independent cross sections. We interpret this 
as being due to the fact that a large spin-dependent cross section implies efficient 
scattering on hydrogen in the Sun, whereas a large spin- independent cross sec- 
tion is necessary for capture in the Earth (which is composed of spinless nuclei). 
The highest fluxes also occur for rather light masses of the CP-odd Higgs boson, 
m H o < 150 GeV. 
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Figure 5: The muon flux obtained from a) only the old population and b) both the old and 
the new population of WIMPs. The points are symbol coded according to whether they 



exceed the claimed direct detection bound a si [16] from the Dam A experiment (squares), 
give direct detection rates within a factor of 10 of those limits (pluses) or give rates more 
than 10 times smaller than these limits (crosses). 



In future experiments, it may be realistic to reach a sensitivity of the order of 
a few events per km 2 per year, covering most of the area in the flux versus mass 
diagram in fig. |5|. It is obvious that the enhancements found in this work may be 
crucial if, for instance, the uncorrected rate would fall just below this value. (It is to 
be reminded that a priori, we have no knowledge whatsoever about which model, if 
any, would represent "the world" in our scans.) 

The correlation between signals from the Earth and the Sun have also been 
investigated. Even without this new population, there are models for which the flux 
from the Earth is higher than that from the Sun, but including this new population, 
we find many more models (at 60-130 GeV) for which it is more advantageous to 
look at the Earth than the Sun. 

We have also made runs with the angular momentum parameter A = 2. The 
results are very similar, and we do not display them here. The main effect is a shift 
of the kinematical cut-off at 150 GeV to around 170 GeV, and an increased muon 
rate by up to 50 % in some cases. 
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8. Conclusions 

We have investigated the contribution from a new population of WIMPs coming 
from WIMPs that have scattered in the outskirts of the Sun and (due to small 
perturbations by the other planets) are confined to nearly radial orbits reaching out 
to the Earth. 

We find that even for the conservative case when A = 1 (i.e. the z component of 
the angular momentum is conserved), we can get enhancements of up to a factor of 
100 when the neutralino mass is less than around 150 GeV. 

For the mass range of 60-130 GeV some models can already be ruled out on the 
basis of existing data. The next generation of detectors should be able to probe a 
much larger region of parameter space in this mass range than would otherwise be 
possible. 

As a result of this new WIMP population, for certain models with WIMPs in 
the mass range 60-130 GeV, there can be an even greater advantage than before in 
looking for the indirect signal from the Earth rather than from the Sun. 

Whether enhanced annihilations in the Earth for larger WIMP masses might be 
possible depends on knowing the precise details of the Solar System WIMP velocity 
distribution, which will await further numerical work on the evolution of Solar System 
WIMPs. 
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